#include <algorithm>
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <unistd.h>

typedef std::vector< std::pair<int, std::vector<unsigned int> > > Matrix;

std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems) {
    std::stringstream ss(s);
    std::string item;
    while(std::getline(ss, item, delim)) {
        elems.push_back(item);
    }
    return elems;
}

std::vector<std::string> split(const std::string &s, char delim) {
    std::vector<std::string> elems;
    return split(s, delim, elems);
}

int str2int(std::string str) {
    return atoi(str.c_str());
}

float str2float(std::string str) {
    return atof(str.c_str());
}

int main( int argc, char *argv[] ) {


    int opt = 0;
    int genomes_num = -1;
    const char *matrixFileName = NULL;    /* параметр -m */
    const char *vectorFileName = NULL;    /* параметр -v */
    static const char *optString = "g:m:v:";
    
    // read options
    opt = getopt( argc, argv, optString );
    while( opt != -1 ) {
        switch( opt ) {
            case 'g':
                genomes_num = atoi(optarg);
                break;
            case 'm':
                matrixFileName = optarg;
                break;
            case 'v':
                vectorFileName = optarg;
                break;
            default:
                /* сюда на самом деле попасть невозможно. */
				break;
        }
        opt = getopt( argc, argv, optString );
    }    

    std::ifstream matrix_file (matrixFileName);
    if (not matrix_file.is_open()) {
        std::cerr << "Unable to open matrix file" << std::endl;
        return 1;
    }
    matrix_file.close();

    std::ifstream vectors_file (vectorFileName);
    if (not vectors_file.is_open()) {
        std::cerr << "Unable to open vectors file " << std::endl;
        return 1;
    }
    
    if (genomes_num <=0) {
        std::cerr << "Genomes_num is incorrect: " << genomes_num << std::endl;
	return 1;
    }
    
    std::vector<std::string> v_strs;
    std::vector<std::string> v_gene_ids_str;
    std::vector<int> v_gene_ids;
    std::vector<std::string> m_strs;
    std::vector<std::string> m_gene_ids_str;
    std::vector<int> m_gene_ids;
    std::map<int, float> m_gene_ids_weights_map;
    std::map< std::pair<int, int>, double> probability_cache;
    std::vector<std::string>::iterator id_weight_iter;
    std::vector<std::string> id_weight;
    std::vector<int>::iterator set_it;
    
    int n1, n2, n12, genom_id;
    float nf2, nf12, weight;
    double probability;
    std::pair<int,int> cache_key;
    std::string vline;
    std::string mline;

    while (vectors_file.good()) {
        getline(vectors_file, vline);
        if (vline.size() == 0) { break;	}
        // parse line from vector file
        v_strs = split(vline, '\t');
        v_gene_ids_str = split(v_strs[1], ',');
        v_gene_ids.resize(v_gene_ids_str.size());
        std::transform(v_gene_ids_str.begin(), v_gene_ids_str.end(), v_gene_ids.begin(), str2int);
        std::sort(v_gene_ids.begin(), v_gene_ids.end());
        n1 = v_gene_ids.size();

        probability_cache.clear();
        matrix_file.open(matrixFileName);
        while (matrix_file.good()) {
//            std::cout << "read matrix" << std::endl;
            getline(matrix_file, mline);
            if (mline.size() == 0) { break;	}
            // parse line from matrix file
            m_strs = split(mline, '\t');
            m_gene_ids_str = split(m_strs[1], ',');
            m_gene_ids.clear();
            m_gene_ids_weights_map.clear();
            for (id_weight_iter=m_gene_ids_str.begin(); id_weight_iter!=m_gene_ids_str.end(); ++id_weight_iter) {
            	id_weight = split(*id_weight_iter, ':');
//                std::cout << weight << std::endl;
                weight = str2float(id_weight[1]);
                genom_id = str2int(id_weight[0]);
                if (weight > 0.0) {
                    m_gene_ids.push_back(genom_id);
                    m_gene_ids_weights_map[genom_id] = weight;
                }
            }
            std::sort(m_gene_ids.begin(), m_gene_ids.end());
            // intersect genom_ids
            std::vector<int> gene_ids_set;
            std::insert_iterator<std::vector<int> > gene_ids_set_ins(gene_ids_set, gene_ids_set.begin());
            std::set_intersection(m_gene_ids.begin(), m_gene_ids.end(), v_gene_ids.begin(), v_gene_ids.end(), gene_ids_set_ins);
            nf12 = 0.5;
            for (set_it=gene_ids_set.begin(); set_it!=gene_ids_set.end(); set_it++){
                nf12 += m_gene_ids_weights_map[*set_it];
            }
            n12 = int(nf12);
            nf2 = 0.5;
            for (set_it=m_gene_ids.begin(); set_it!=m_gene_ids.end(); set_it++){
                nf2 += m_gene_ids_weights_map[*set_it];
            }
            n2 = int(nf2);
            // print output
            if (n12 > float(n1*n2)/genomes_num) { // больше среднего
//            if (true) {
                cache_key = std::make_pair(n2, n12);
                if (probability_cache.count(cache_key)>0) {
                    probability = probability_cache[cache_key];
                }
                else {
                    probability = theor_probability_gt(genomes_num, n1, n2, n12);
                    probability_cache[cache_key] = probability;
                }
                if (probability <= 0.01) {
         	        std::cout   << v_strs[0] << "\t" << m_strs[0] << "\t" \
                        << probability << "\t" \
                        << n1 << "," << n2 << "," << n12 \
                        << std::endl;
                }
            }
        }
        matrix_file.close();
    }
    vectors_file.close();
    return 0;

}
